###############################################
###                                         ###
###              Read the Data              ###
###                                         ###
###############################################
#### Clean Up
rm(list=ls())
setwd("~/Dropbox/Shanshan Lian/ISA_MIDWEST 2020 (Co-Authored Paper)/Submission to JHR/RR JHR/Replication File_Update")
library("geepack")
library("dplyr")
library("foreign")
Data_All <- read.csv("Data_All.csv")
Data_model1 <- read.csv("Data_model1.csv")
Data_model2 <- read.csv("Data_model2.csv")
Data_model3 <- read.csv("Data_model3.csv")
Data_model4 <- read.csv("Data_model4.csv")
Data_icews <- read.csv("Data_ICEWS.csv")


##########################################################################
###                                                                    ###
###                  Tables and Figures in the Manuscript              ###
###                                                                    ###
##########################################################################
###############################################
###                                         ###
###                 Table 1                 ###
###                                         ###
###############################################

### Model 1
model1 <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2xcs_ccsi, 2)+ 
                  p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                  data=Data_model1)
summary(model1)            

### Model 2
model2 <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2cseeorgs, 2)+
                 p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                 data=Data_model2)
summary(model2)

### Model 3
model3 <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2csreprss, 2)+
                 p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                 data=Data_model3)
summary(model3)

### Model 4
model4 <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2csprtcpt, 2)+
                 p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                 data=Data_model4)
summary(model4)

###############################################
###                                         ###
###              Figures 4 & 5              ###
###                                         ###
###############################################

library("ggplot2")
library("jtools")
library("gridExtra")
### Figures 4
plot_model1<- effect_plot(model1, pred=v2xcs_ccsi, data=Data_model1,interval=TRUE,
                          x.label="Core Civil Society Index from V-Dem 
            (Scale Reversed so Higher Values = More Repressive Environment)",
                          y.label = "Predicted NGO-to-Government 
            Conflict-to-Cooperation Interactions")
plot_model1

### Figures 5
plot_model2 <- effect_plot(model2, pred=v2cseeorgs, data=Data_model2,interval=TRUE,
                         x.label="CSO Entry and Exit from V-Dem 
                         (Scale Reversed so Higher Values = More Repressive Environment)",
                         y.label = "Predicted NGO-to-Government 
            Conflict-to-Cooperation Interactions")


plot_model3 <- effect_plot(model3, pred=v2csreprss, data=Data_model3,interval=TRUE,
                         x.label="CSO Repression from V-Dem 
                        (Scale Reversed so Higher Values = More Repressive Environment)",
                         y.label = "Predicted NGO-to-Government 
            Conflict-to-Cooperation Interactions")

plot_model4 <- effect_plot(model4, pred=v2csprtcpt, data=Data_model4,interval=TRUE,
                         x.label="CSO Participatory Environment from V-Dem 
            (Scale Reversed so Higher Values = More Repressive Environment)",
                         y.label = "Predicted NGO-to-Government 
            Conflict-to-Cooperation Interactions")

grid.arrange(plot_model2, plot_model3,plot_model4, ncol=2)

##########################################################################
###                                                                    ###
###               Tables and Figures in the Online Appendix            ###
###                                                                    ###
##########################################################################
###############################################
###                                         ###
###                 Table 1                 ###
###                                         ###
###############################################
Data_All %>% 
  dplyr::select(F.meanNGOtoGovtIntensityICEWS, meanNGOtoGovtIntensityICEWS,
                v2xcs_ccsi, v2cseeorgs, v2csreprss, v2csprtcpt, lnicews_repression_intensity,
                p_polity2, lnAidAllNGOs, lnNonNGOAid, KOFPoGIdf, lnhrcountfilled) %>% 
  data.frame() %>% 
  stargazer::stargazer(header = FALSE,
                       title = "Summary of the Numeric Variables",
                       type = "text")

###############################################
###                                         ###
###                 Figure 1                ###
###                                         ###
###############################################
hist(Data_All$F.meanNGOtoGovtIntensityICEWS, main=" ",
     xlab="Mean Intensity of NGO-to-Government Conflictual-to-Cooperative Interactions (t+1) ",
     xlim=c(-7,7), ylim=c(0,3500),
     col="lightpink1")

###############################################
###                                         ###
###              Figure 2 & 3               ###
###                                         ###
###############################################

### Figure 2
# model 1
loess_model1 <- loess(F.meanNGOtoGovtIntensityICEWS~v2xcs_ccsi, data=Data_model1)
smooth_model1 <- predict(loess_model1, se = TRUE)

plot(F.meanNGOtoGovtIntensityICEWS~v2xcs_ccsi,data=Data_model1, 
     ylim=c(-0.8,0.4),pch=1, cex=0.5, col="grey",
     xlab="Core Civil Society Index from V-Dem (Scale Reversed so Higher Values = More Repressive Environment)", ylab="NGO-to-Government Conflict-to-Cooperation Interactions (t+1)")
a <- order(Data_model1$v2xcs_ccsi)
lines(Data_model1$v2xcs_ccsi[a], smooth_model1$fit[a], col="red", lwd=3, lty=1)
lines(Data_model1$v2xcs_ccsi[a], smooth_model1$fit[a] - qt(0.025, smooth_model1$df)*smooth_model1$se[a],lty=1,lwd=1, col="red")
lines(Data_model1$v2xcs_ccsi[a], smooth_model1$fit[a] - qt(0.975, smooth_model1$df)*smooth_model1$se[a],lty=1, lwd=1,col="red")


### Figure 3
par(mfrow=c(2,2))

# model 2
loess_model2 <- loess(F.meanNGOtoGovtIntensityICEWS~v2cseeorgs, data=Data_model2)
smooth_model2 <- predict(loess_model2, se = TRUE)
plot(F.meanNGOtoGovtIntensityICEWS~v2cseeorgs,data=Data_model2, 
     ylim=c(-0.8,0.4),pch=1, cex=0.5, col="grey",
     xlab="CSO Entry and Exit from V-Dem (Scale Reversed so Higher Values = More Repressive Environment)", ylab="NGO-to-Government Conflict-to-Cooperation Interactions (t+1)")
a <- order(Data_model2$v2cseeorgs)
lines(Data_model2$v2cseeorgs[a], smooth_model2$fit[a], col="red", lwd=3, lty=1)
lines(Data_model2$v2cseeorgs[a], smooth_model2$fit[a] - qt(0.025, smooth_model2$df)*smooth_model2$se[a],lty=1,lwd=1, col="red")
lines(Data_model2$v2cseeorgs[a], smooth_model2$fit[a] - qt(0.975, smooth_model2$df)*smooth_model2$se[a],lty=1, lwd=1,col="red")

# model 3
loess_model3 <- loess(F.meanNGOtoGovtIntensityICEWS~v2csreprss, data=Data_model3)
smooth_model3 <- predict(loess_model3, se = TRUE)
plot(F.meanNGOtoGovtIntensityICEWS~v2csreprss,data=Data_model3, 
     ylim=c(-0.8,0.4),pch=1, cex=0.5, col="grey",
     xlab="CSO Repression from V-Dem (Scale Reversed so Higher Values = More Repressive Environment)", ylab="NGO-to-Government Conflict-to-Cooperation Interactions (t+1)")
a <- order(Data_model3$v2csreprss)
lines(Data_model3$v2csreprss[a], smooth_model3$fit[a], col="red", lwd=3, lty=1)
lines(Data_model3$v2csreprss[a], smooth_model3$fit[a] - qt(0.025, smooth_model3$df)*smooth_model3$se[a],lty=1,lwd=1, col="red")
lines(Data_model3$v2csreprss[a], smooth_model3$fit[a] - qt(0.975, smooth_model3$df)*smooth_model3$se[a],lty=1, lwd=1,col="red")

# model 4
loess_model4 <- loess(F.meanNGOtoGovtIntensityICEWS~v2csprtcpt, data=Data_model4)
smooth_model4 <- predict(loess_model4, se = TRUE)
plot(F.meanNGOtoGovtIntensityICEWS~v2csprtcpt,data=Data_model4, 
     ylim=c(-0.8,0.4),pch=1, cex=0.5, col="grey",
     xlab="CSO Participatory Environment from V-Dem (Scale Reversed so Higher Values = More Repressive Environment)", ylab="NGO-to-Government Conflict-to-Cooperation Interactions (t+1)")
a <- order(Data_model4$v2csprtcpt)
lines(Data_model4$v2csprtcpt[a], smooth_model4$fit[a], col="red", lwd=3, lty=1)
lines(Data_model4$v2csprtcpt[a], smooth_model4$fit[a] - qt(0.025, smooth_model4$df)*smooth_model4$se[a],lty=1,lwd=1, col="red")
lines(Data_model4$v2csprtcpt[a], smooth_model4$fit[a] - qt(0.975, smooth_model4$df)*smooth_model4$se[a],lty=1, lwd=1,col="red")

dev.off()

###############################################
###                                         ###
###                 Table 2                 ###
###                                         ###
###############################################
model_icews <- lm(F.meanNGOtoGovtIntensityICEWS~poly(lnicews_repression_intensity, 2)+
                  p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                data=Data_icews)
summary(model_icews)


###############################################
###                                         ###
###                Figure 4                 ###
###                                         ###
###############################################
plot_icews <- effect_plot(model_icews, pred=lnicews_repression_intensity, data=Data_icews,interval=TRUE,
                        x.label="CSO Repression from ICEWS (Higher Values = More Repressive Environment)",
                        y.label = "Predicted NGO-to-Government 
            Conflict-to-Cooperation Interactions")
plot_icews


###############################################
###                                         ###
###                 Table 3                 ###
###                                         ###
###############################################

library("lmtest")
library("sandwich")

coeftest(model1, vcov = vcovHC(model1, type = "HC1"))
coeftest(model2, vcov = vcovHC(model2, type = "HC1"))
coeftest(model3, vcov = vcovHC(model3, type = "HC1"))
coeftest(model4, vcov = vcovHC(model4, type = "HC1"))

###############################################
###                                         ###
###                 Table 4                 ###
###                                         ###
###############################################

### Model 1
model1_lag <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2xcs_ccsi, 2)+  meanNGOtoGovtIntensityICEWS+
                     p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                     data=Data_model1)
summary(model1_lag)

### Model 2
model2_lag <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2cseeorgs, 2)+meanNGOtoGovtIntensityICEWS+
                    p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                    data=Data_model2)
summary(model2_lag)

### Model 3
model3_lag <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2csreprss, 2)+meanNGOtoGovtIntensityICEWS+
                    p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                    data=Data_model3)
summary(model3_lag)

### Model 4
model4_lag <- lm(F.meanNGOtoGovtIntensityICEWS~poly(v2csprtcpt, 2)+meanNGOtoGovtIntensityICEWS+
                    p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                    data=Data_model4)
summary(model4_lag)


###############################################
###                                         ###
###                 Table 5                 ###
###                                         ###
###############################################

### Model 1
Data_model1_gee <- na.omit(Data_model1)
model1_gee <- geeglm(meanNGOtoGovtIntensityICEWS~poly(v2xcs_ccsi, 2)+
                            p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                    data=Data_model1_gee,
                    family=gaussian,
                    id=iso, corstr='ar1')
summary(model1_gee)

### Model 2
Data_model2_gee <- na.omit(Data_model2)
model2_gee <- geeglm(meanNGOtoGovtIntensityICEWS~poly(v2cseeorgs, 2)+
                           p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                   data=Data_model2_gee,
                   family=gaussian,
                   id=iso, corstr='ar1')
summary(model2_gee)

### Model 3
Data_model3_gee <- na.omit(Data_model3)
model3_gee <- geeglm(meanNGOtoGovtIntensityICEWS~poly(v2csreprss, 2)+
                           p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                   data=Data_model3_gee,
                   family=gaussian,
                   id=iso, corstr='ar1')
summary(model3_gee)

### Model 4
Data_model4_gee <- na.omit(Data_model4)
model4_gee <- geeglm(meanNGOtoGovtIntensityICEWS~poly(v2csprtcpt, 2)+
                           p_polity2+lnAidAllNGOs+lnNonNGOAid+OECD+KOFPoGIdf+lnhrcountfilled+locationconflict, 
                   data=Data_model4_gee,
                   family=gaussian,
                   id=iso, corstr='ar1')
summary(model4_gee)























